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Abstract 

In this paper the combination of the Osher approxi- 
mate Riemann solver for the Euler equations and vari- 
ous ENO schemes is discussed for one-dimensional flow, 
The three basic approaches, viz. the ENO scheme us- 
ing primitive variable reconstruction, either with the 
Cauchy-Kowalewski procedure for time integration or 
theTVD Runge-Kutta scheme, and the flux-ENO meth- 
od are tested on different shock tube cases. The shock 
tube cases were chosen to present a serious challenge 
to the ENO schemes in order to test their ability to 
capture flow discontinuities, such as shocks. Also the 
effect of the ordering of the eigenvalues, viz. natural or 
reversed ordering, in the Osher scheme is investigated. 
The ENO schemes are tested up to fifth order accu- 
racy in space and time. The ENO-Osher scheme using 
the Cauchy-Kowalewski procedure for time integration 
is found to be the most accurate and robust compared 
with the other methods and is also computationally effi- 
cient. The tests showed that the ENO schemes perform 
reasonably well, but have problems in cases where tw'o 
discontinuites are close together. In that case there are 
not enough points in the smooth part of the flow to 
create a non-oscillatory interpolation. 

1. Introduction 

The development of high order accurate, non-oscillatory 

shock capturing schemes currently is an area of active 

interest. High order accuracy is important for more 

complicated unsteady inviscid problems and for direct 
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simulation of compressible flows. It is fairly straightfor- 
ward to incorporate high order accuracy in non-conser- 
vative finite difference methods, however, shock cap- 
turing will not be possible. Finite volume methods 
and conservative finite difference methods, which re- 
tain this property, are unfortunately limited to first or 
second order accuracy in most cases. An important rea- 
son for this limitation in accuracy is the use of Total 
Variation Diminishing (TVD) methods to obtain non- 
oscillatory solutions. TVD methods are limited to first 
order accuracy in more than one dimension, Goodman 
and Leveque 1 * 3 , and even in one dimension they reduce 
to first order accuracy at non-sonic local extrema, Os- 
her and Chakravarthy 11 * * . 

Harten, Osher, Shu et al . 4,5,0 * 7,14 * 15 developed in recent 
years the so-called Essentially Non-Oscillatory (ENO) 
schemes, which do not have this limitation and have 
uniform high order accuracy outside discontinuities. In 
this paper the combination of the Osher approximate 
Riemann solver and the various ENO schemes for the 
solution of the Euler equations will be discussed. The 
discussion will be limited to one-dimensional problems, 
but important information for the development of multi- 
dimensional ENO methods can be obtained from the 
solution of several shock tube problems. 

The main feature of ENO schemes is that they use an 
adaptive stencil. At each grid point a searching al- 
gorithm determines which part of the flow surround- 
ing that grid point is the smoothest. This stencil is 
then used to construct a higher order accurate, conser- 
vative interpolation to determine the variables at the 
cell faces. This interpolation process can be applied 
to the conservative variables, characteristic variables 
or the fluxes, either defined as cell averaged or point 
values. The ENO scheme tries to minimize numerical 
oscillations around discontinuities by using predomi- 
nantly data from the smooth parts of the flow field. 
Due to the constant stencil switching the ENO scheme 
is highly non-linear and only limited theoretical results 
are available, Harten and Osher 5 * * * * * and Harten et al. 4 . 
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An important subject in the ENO schemes is the nu- 
merical approximation of the fluxes at the cell faces. 
The Roe approximate Riemann solver is frequently used, 
while the Lax-Friedrichs flux splitting is preferred when 
the differentiability of the flux is important. The Roe 
scheme, however, requires a modification, because it 
does not satisfy the entropy condition and allows steady 
expansion shocks, whereas the Lax-Friedrichs scheme is 
very dissipative. When using an exact Riemann solver 
instead of the Roe approximation Harten et al. 6 ob- 
tained better results, while also the results of Shu and 
Osher 14,15 improved with a better flux approximation. 

The Osher scheme, which is the most accurate approx- 
imate Riemann solver available, has some nice proper- 
ties, which might improve the accuracy of ENO schemes. 
Until now the combination of the OsheT and ENO scheme 
has not been investigated. Compared with the Roe 
approximate Riemann solver the Osher scheme is less 
popular due to the higher cost of computing t lie fluxes 
and the fact that it is more complicated than the Roe 
scheme. When using a higher order ENO scheme the 
cost of the Riemann solver, however, becomes less dom- 
inant, because the ENO interpolations themselves are 
quite expensive, especially for multi-dimensional prob- 
lems. In addition a careful implementation of the Osher 
scheme does not have to be too expensive. The use of 
the Osher flux splitting has some additional benefits. It 
has a very low numerical dissipation in boundary lay- 
ers compared to other Euler schemes, see for instance 
Koren 8 . Although beyond the scope of this papei this 
is an important consideration in the approximation of 
the inviscid contribution when solving the compress- 
ible Navier-Stokes equations for flows with both strong 
shocks and boundary layers by means of direct simu- 
lations. The possible application of ENO schemes in 
direct numerical simulations of compressible flow is the 
reason that the ENO schemes are tested up to fifth or- 
der accuracy. 

As already mentioned there ate various ways to con- 
struct ENO schemes. A good comparison of the dif- 
ferent approaches unfortunately is lacking and will be 
part of this paper, with special emphasis on the combi- 
nation of the Osher and ENO schemes. Several difficult 
shock tube cases are computed, which present a serious 
challenge to ENO schemes. In addition the effect of the 
ordering of the eigenvalues in tlie Osher scheme will be 
investigated, which is one of parameters which might 
effect the accuracy of the schemes. This comparison 
gives valuable information which of the approaches is 
the most successful. In the next section first a brief dis- 
cussion of the Osher scheme for the Euler equat ions will 
be given. Next the construction of the various Osher- 


ENO schemes will be discussed and results of several 
shock tube tests will be given at the end of the paper. 

2. Osher Scheme 

The Osher-ENO scheme is used to solve the Euler equa- 
tions of gas dynamics, which are defined as: 

u t -F f (u) x — 0 i > 0, -oo < x < oc 

(2-1) 

with initial data: u(z,0) — Uo- The vectors u and 
f are defined as n — (p, pu,e) f and f =• (pu, pu 2 , (e -f 
p)u)F Here p, it, e and p represent the density, velocity, 
total energy and pressure respectively. The Jacobian 
matrix df has real eigenvalues A T , (A 1 = u — c, v, u 4- c), 
with c the speed of sound, and has a complete set of 
eigenvectors r*. 

The system of equations (2.1) is completed with an 
equation of state for an ideal gas, p =• pRT ) with R 
the gas constant, and the relation between total energy 
and temperature, T = e/(pc v ) — \u 2 , where c v repre- 
sents the specific heat at constant volume. 

Let uj 1 = u(xi,t n ), x i ~ jh\t n = nr , with grid spacing 
h and time step r, then the conservative first order 
approximation to the Euler equations can be written 
as: 

< +l =u?-A(h lH - hi _p (2.2) 

with A = r/h. The flux h at a cell face with index i> + \ 
in the Osher scheme now is defined as: 

h t+ i = + fi+l - J r (dr + - at-) du ) (2.3) 

with r t a path in phase space. The matrices df± are 
defined as: df± — 5A ± 5“ 1 , with S -1 the matrix with 
right eigenvectors r* as columns, and 
A + = diag(max(0, A*)) ; — diag(min(0, A 1 )). Here 

diag is short for diagonal matrix. In the construction 
of the flux- ENO scheme we will also use the expression 
for the flux differences: 

h i+ i —h t _i — f df + du-\- f df ~dv (2.4) 

V P t _ I V F , 

By choosing a specific path in phase space, namely 
- — r^s), Osher was able to derive explicit expres- 
sions for the integrals in equations (2.3-4) by separating 
the integrals along T l into three parts. Along each path 
one eigenvalue A* is constant and has a related set 
of Riemann invariants, A 1,3 : ^ = {u ± ^7) and 

A 2 : ip — {u,p}. These Riemann invariants ^ are used 
to determine the intermediate states by linking them 
to Uj and Uj+i, which are the end points of the path 
F,. The order in which the eigenvalues are used in the 
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computation of the path integrals along T, is impor- 
tant. The order (u - c, u, u + c) is called the natural 
order, whereas (u-f c, u } u — c) is the reversed order. Os- 
her and Solomon 9 were able to prove that the reversed 
order admits steady shocks in at most two points and 
satisfies the entropy condition, which guarantees con- 
vergence almost everywhere to the solution of the Euler 
equations in the limit h — ♦ 0, r — ♦ 0. 

The integrals along each path can be computed using 
the fact that the eigenvalues are either genuinely non- 
linear, (V U A* * r* = 1), or linearly degenerate, (V U A* • 
r* = 0). In the first case A can have at most one zero 
on each path and in the second case it has no zeros. If 
there is a zero on a path, then the path integral has 
to be modified to treat the sonic point. The careful 
treatment of sonic points eliminates the possibility of 
expansion shocks which appear in the Roe method and 
also gives a differentiable flux vector. 

The Osher scheme uses simple waves to solve the Rie- 
mann problem, which become multivalued in a shock. 
For an analysis of this phenomenon see Van Leer 16 
and Dubois 2 . More details about the implementation 
of the Osher scheme and explicit expressions for the 
flux integrals can be found in Osher and Solomon 9 , 
Chakravarthy and Osher 1 and Rai and Chakravarthy 13 . 

3. ENQ Sc hemes 

ENO schemes overcome the limitations of TVD schemes 
by relaxing the requirement of total variation non-increa- 
sing. They are conservative, essentially non-oscillatory 
and give uniform accuracy in smooth regions, without 
the degradation of accuracy at non-sonic local extrema 
as observed with TVD methods. There are several 
approaches possible when constructing ENO schemes. 
Harten, Osher, Engquist and Chakravarthy 6 use the 
ENO scheme to construct a higher ordeT solution to 
the cell-aveTage of equation (2.1) using a sliding aver- 
age, defined as: 

1 f * 

u(x >0 = W u[x + y,i)dy (3.1) 

Integrating equation (2.1) over the domain x t+ i] 

gives: 

^u(x,<)+i[f(u(® + |,t))-f(u(x-^,t))] = 0 (3.2) 

with grid spacing h. When considering this equation at 
the point x — aq, the cell center, this gives the method 
of lines formulation for the cell average u (x l} t). Inte- 
grating in time finally gives the cell averaged equation: 

u(xi,< + r) = u(xi,<) - A[h(x, + 1 , t; u) - h(x ; _ . , t; u)] 

3 (3.3) 


with A — r/h and 


h(x, <; u) = - f f(u (x,t + ii))dr] (3.4) 

T Jo 

The numerical flux in the ENO scheme from Harten et 
al. now is constructed such that it approximates the 
exact flux up to 0(h r ): 

u ? +I -fE fc (r)-u"]. = 0(A T ) (3.5) 

with E h {r) the numerical solution operator. The ENO 
scheme of Harten et al. therefore gives an r-th order 
accurate approximation to the cell averages. The most 
important ingredient of their ENO method is the recon- 
struction of the point values n(x) from the cell averaged 
values Ui. These point values are necessary to compute 
the flux i i at the cell faces. This is done with a 

-j 

reconstruction method, discussed in the next section, 
such that the piecewise polynomial R(x, u) is conser- 
vative, viz. R(x{, u) = u,, essentially non-oscillatory, 
TV(R{ • ; u)) < TV (u) + 0(/i r ), and gives at all points 
in a neighborhood around x, an r-th order approxima- 
tion to u, when u is smooth. The time integral, equa- 
tion(3.4), then is computed using the Cauchy-Kowalew- 
ski procedure for hyperbolic systems and is discussed 
in the next section. This method has as benefit that it 
couples the spatial and temporal discretization. 

In an alternative approach Shu and Osher 14,15 con- 
structed the numerical solution directly from equation 
(3.2), with x — Xi, using a correction to the fluxes, 
obtained by a Taylor series expansion around the cell 
center x^. 




+ ^2 a 2kh 2k 


fe =1 


/<9 2fc h\ 
\d x& 




+ 0(/i 2m+1 ) 


(3 - 6) 

Here h t+ i is the flux obtained with any first order con- 
servative scheme. The first few coefficients are: a2 — 
-^,a4 = 5^501 ’*’• This procedure makes the conser- 
vative formulation higher order accurate. This method 
is essentially a finite difference method and not a finite 
volume method, because u(x X) t) is now considered as a 
point value in the Taylor series expansion instead of a 
cell average. It is a generalization of the high resolution 
TVD schemes proposed by Osher and Chakravarthy 12 , 


In their second paper Shu and Osher 15 demonstrated 
that it is possible to construct a higher order method by 
reconstructing the fluxes directly from the fluxes com- 
puted from the cell averaged values without using this 
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formula. In this case the variable u should be consid- 
ered again as a cell averaged variable. This approach 
is used in section 3.2 for the flux-ENO method. 

The time integration is performed using a new Runge- 
Kutta scheme which does not increase the total vari- 
ation in time. This approach has as benefit that the 
spatial and temporal discretization are decoupled. This 
results in a simpler coding, but increased cost because 
the reconstruction, which is the most expensive part of 
the ENO schemes, has to be done at each stage of the 
Runge-Kutta scheme. 

A hybrid approach is also possible, namely using the 
ENO scheme from Harten et al., but with the TVD 
Runge-Kutta method from Shu and Osher instead of 
the Cauchy- Kowalewski procedure for time integration. 
This method saves the cumbersome derivation of the 
higher derivatives of the solution but still does the ENO 
reconstruction on the primitive variables instead of the 
fluxes, which is generally more accurate. Results of this 
approach will be presented at the end of the paper. 
The main benefit of the original point ENO method 
of Shu and Osher is that it is easy to extend to more 
dimensions using dimension splitting. Although in this 
way a higher order method, is obtained in practice it 
has a strong bias to the principle directions of the grid, 
which seriously degrades the accuracy. The approach 
of Harten et al. will give uniform higher order accuracy 
with a minimal grid dependence, but at significant cost 
and complexity, see Harten and Chakravarthy 7 . 

3.1 ENO- Osher scheme based on reconstruction from 
cell averaged variables 

The implementation of the Osher scheme in the ENO 
method from Harten et al. 6 , which uses a reconstruc- 
tion from the cell averaged variables is straightforward. 
Only the Riemann solver has to be replaced with the 
Osher approximate Riemann solver, discussed in sec- 
tion 2. For completeness, however, a summary of this 
algorithm will be given. More details can be found in 
Harten et al. 0 . The first step in the ENO reconstruction 
is the determination of the primitive function JJ(x t+ i) 
from the cell averaged variables: 

i 

u ( x »+|) = X ^ x '+5 “ ( 3 . 1 . 1 ) 

j=o 

A higher order polynomial representation of U in each 
cell is now constructed by determining the divided dif- 
ferences used in the Newton interpolation method using 
the following recursive algorithm: Start with k min — 
i — | and k max = i -f and the divided differences: 

H[x,_i] = U(*,_i); H[,T;_ i , x i + i] (3.1.2) 


If the divided difference x i+ j] is larger 

than H[a^_ i , x i+%] choose 

i , z t +|] and k mai = otherwise 

H[x t _ $ , x i _i, ^ ] is accepted and k min = x t _s. This 
process is repeated till the required order of the interpo- 
lation is obtained and applied to each component of U 
independently. It will give the interpolation of U and 
the stencil [k m i n ,k max ]. Here the divided differences 
are defined as: 

' * * > x i+k + i] — , * ' ' > + 

+ , * ’ ■ , £t + fc_i])/(£t + fc + i — X i + ± ) 

(3.1.3) 

After the determination of the coefficients of the New- 
ton polynomial it is straightforward to obtain the point 
values by differentiation of the Newton interpolation 
polynomial, because u = ^U. This process is greatly 
simplified using the algorithm discussed in the Appendix 
in Harten et al. 6 , which transforms the Newton poly- 
nomial representation into a series expansion around 
.r t : 

r - 1 

J?(x;u) = y^b, [t (x - Xj) k (3.1.4) 

k=0 

This process gives a representation of the solution in 
each cell and can be used to determine the values of 
u at the cell faces The values at the left and 

right side of the cell are now used in the Osher ap- 
proximate Riemann solver, which gives the fluxes $ t ±i- 
This reconstruction process can be applied directly to 
the conservative variables, but in order to minimize the 
interactions between the different equations it is prefer- 
able first to transform u to characteristic variables and 
transform back after the reconstruction. 

The method of lines approach can now be followed, 
solving equation (3.2), with the TVD Runge-Kutta met- 
hod, or the Cauchy-Kowalewski procedure can be used 
to solve equation (3.3). The TVD Runge-Kutta method 
is an r-th order, r stage Runge-Kutta method defined 
as: 


» — 1 


U (,) = ^ + Aj trL(u (fc) )] , 

i ~ 1 , • • - r 

k-0 

(3.1.5) 

with 


u (0) _ u » u (r) _ u n + l 

(3.1.6) 


Here L represents the spatial discretization operator. 
For the coefficients a, (3 see Shu and Osher 14 . The 
TVD Runge-Kutta method has as main benefit that it 
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does not increase the total variation during the time 
integration step. 

The alternative is to use equation (3.3) and compute 
the integral, equation (3.4). This integral is discretized 
by means of a Gaussian quadrature formula: 



K 

t+Tj))dr) - (u(x i+ . 

k= 0 


Ar))+0(r') 


(3-1.7) 

The values of u(x 1 + i , Pk^) are now determined using 
the Cauchy-Kowalewski procedure, which is essentially 
a Taylor series expansion around the point (x,,/ 71 ): 


r — 1 I 


u(x,t) = 

1=0 k = o 


d*u(Xi, t n ) (x — Xi) k t l ~ k 

dx k dt l - k hi (l-k)\ 


(3.1.8) 


The derivatives can be obtained by differen- 

tiation of the original differential equation (2.1) and 
__ ■ * £1 1 
using the coefficients which are equal to 


<9u(x t ,/ n ) 

dx l 


= b. 


0 < / < r- 1 (3.1.9) 


3.2 Osher fluz-ENO meth od 

The Osher flux-ENO method is constructed using the 
flux difference relation for the Osher scheme, equation 
(2.4). The ENO reconstruction is applied directly to 
the fluxes in a conservative way, as proposed by Shu and 
Osher 15 . The reconstruction is, however, more compli- 
cated due to the path integrals in the Osher flux, equa- 
tion (2.3). Define the flux f in a cell centered around 
x t as f = f g(£) then f satisfies the relation 

f(uO = F(* j+ j)-F (*i. 4 ))/fc (3-2.1) 

with the primitive function F defined as: 

F(x) ~ / g(0^- The conservative flux difference for 
a cell with index i then can be determined directly from 
the primitive function F, which is determined with the 
ENO reconstruction technique. The flux-ENO recon- 
struction is now applied directly to the positive and 
negative fluxes in the Osher scheme, defined as: 

dff +1 = J d( t du (3.2.2) 

These relations, cannot be used directly because the 
primitive functions F^ are not known, nor the fluxes f ± 
themselves. Their explicit form, however, is not needed, 


only their divided differences which can be linked to 
by means of the following relations: 

df + + . =f, + - f- + _, 

= F+[x i+ , ) x i _,]-F+[x i _,,x i _ f ] (3.2.3) 

with an equivalent relation for df “ , : 

c+i-fr 

F~[x t+ |,x i+ i] - F-[x, + ,,x,_i] (3.2.4) 

(*i+j 

Here F* [xo, ’ ■ * » £*] represent the k-th divided differ- 
ences of F^. These relations automatically introduce 
upwinding for the positive and negative wave direc- 
tions. Higher order approximations can be obtained 
by extending the divided difference tables of F^, us- 
ing the following recursive algorithm for the positive 
flux, which is slightly different from the one discussed 
in section 3.1. 

Start with k mm — i. If the divided difference 
n + [***-7 x ] is larger than 

2 "miu ' 2 

’ ' ' • choose 

Fje * * *» and hrain = ^rmn 1 

irj i u 3 min' 3 

else F. + [x (i_i) i , • • • , x .<i-u 1 1 is accepted and 

* inin *" 2 ? nin "*“* s' 

remains unchanged. This process is done for all the 
components k of the flux vector f + during each in- 
terpolation step l and continued up to the required 
order. The following Newton polynomial now can be 
constructed: 

K min 

Pk~'\*) + 4° n ( x - x -p ( 3 - 2 - 5 ) 

where cj^ is the divided difference accepted in the /- th 
interpolation step. A similar algorithm is used for the 
negative fluxes f~ with k m j n = i, replaced by k m?n = 
i+ 1. Here k m i n is the left most stencil of the grid cen- 
tered around point i\ whereas represents the /-th 
order Newton interpolation polynomial. The index k 
refers to the application of the ENO reconstruction to 
each component of the flux vectors f ^ independently. 
This algorithm automatically chooses the smoothest 
possible interpolation stencil: {k^ n , * ■ • , k^ m -f / — 1}, 
independently for each component of the positive and 
negative fluxes. By continuously comparing the divided 
differences, obtained by adding one point to the left 
and one to the right, it is decided which one gives the 


H'\*) = 
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smoothest interpolation. Adding the positive and nega- 
tive fluxes, and using the relation for the flux difference 
for the Osher scheme, equation (2.4), now gives the flux 
difference for the cell with index i : 

d . 

4 - p+ _ 

+i dx 

(3-2.6) 

The time integration is accomplished using the TVD 
Runge-Kutta scheme of Shu and Osher 14 , discussed in 
the previous section. 

The use of the relations 




- h. i = — P 


dx 


df , + +i - f '+i “ h . + -i 


(3.2.6) 


d(- + ,=h l+ >-f t 

to obtain the positive and negative fluxes, as suggested 
by Shu and Osher 14 , without presenting an application 
in their paper, does not give a higher order scheme. 
The reason is that this relation is based on the average 
of a forward and backward flux, which is only a second 
order accurate approximation to h t + i . In order to ob- 
tain higher order accuracy this average will have to be 
replaced by a higher order interpolation, which would 
require an additional reconstruction step. 

4. Discussion and Results 

To test the combination of the different ENO schemes 
and the Osher scheme several shock tube calculations 
were performed. The cases were chosen to test the 
schemes for the various types of discontinuities which 
exist as solutions of the Euler equations, viz. shocks, 
contact discontinuities and expansion waves. The ini- 
tial conditions were kindly provided by F. Coquel and 
are summarized in Table 1. They were designed to 
be severe tests for shock capturing schemes and are 
helpful in detecting possible flaws in the different ENO 
schemes. For all cases the three possible ENO schemes, 
viz. ENO with primitive variable reconstruction, with 
either the Cauchy-Kowalewski procedure for time inte- 
gration (ENO-CK) or the TVD Runge-Kutta method 
(ENO-RK), or the reconstruction applied directly to 
the fluxes (ENO-FL), were tested. The number of grid 
points was 162 and maximum CFL number .8. In some 
cases the CFL number had to be reduced due to nu- 
merical instabilities, namely case B with ENO-RK and 
case E with ENO-CK required a CFL number .1 for 
4th and 5th order accuracy, whereas for all schemes 
the maximum CFL number had to be reduced to .6 for 
5th order accuracy. All methods were tested with spa- 
tial accuracy ranging from first order up to fifth order. 


The time accuracy was always equal to the spatial accu- 
racy for ENO-CK. For ENO-RK and ENO-FL the time 
accuracy was limited to third order. The TVD Runge- 
Kutta scheme becomes rather awkward for fourth and 
higher order accuracy. It has a strong CFL limitation, 
requires a large amount 
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[N/m 2 

[m/sec] 

[Al 

[N/m 2 ] 

[m/sec] 

[A] 

A 

15000 

0 

1378 

98400 

0 

4390 

B 

988000 

0 

2438 

9930 

0 

2452 

C 

10000 

0 

2627 

100000 

0 

272 

E 

573 

2200 

199 

22300 

0 

546 


Table 1. Initial conditions shock tube cases 


of temporary storage and is not practically useful. In 
addition to testing the various ENO schemes also the 
effect of the ordering of the eigenvalues in the Osher 
scheme has been investigated. For code validation pur- 
poses several of the cases presented by Harten et al. 6 , 
such as Sod's problem were computed and compared 
well with their results. 

The first case, labeled A in Table 1., consists of a left 
moving shock, followed by a contact discontinuity and 
a right moving expansion wave. In all the plots the 
continuous line is the exact solution while the lines with 
dots and triangles represent the numerical results. The 
vertical axis of one curve was shifted upwards to make 
the differences between the various methods more clear. 
All plots, except Fig. 9, show the results using the 
reversed ordering of eigenvalues in the Osher scheme. 
In Fig. 1-3 results of the various methods for 2nd and 
5th order accuracy are presented for the density pro- 
file. The density is the critical variable in this problem, 
because of the contact discontinuity. This means that 
points in the region between shock and contact discon- 
tinuity have a discontinuity in both directions, which 
is the worst possible case for an ENO scheme, because 
in both directions the reconstruction will always have 
to cross a discontinuity for third or higher order re- 
constructions. Due to the extremely low dissipation of 
the higher order schemes these numerical oscillations, 
created in the initial stages when there are not enough 
points in the region between the two discontinuities to 
build a higher order non-oscillatory reconstruction, will 
remain in the solution. It is clear that ENO-CK is su- 
perior in this region, while ENO-RK still gives a rea- 
sonable result, although the method is slightly more 
dissipative than ENO-CK and has a small overshoot. 
The ENO-FL reconstruction, however, has strong oscil- 
lations for third and higher order in the region between 
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the shock and contact discontinuity. Pressure and ve- 
locity are monotone for all schemes and different or- 
ders of accuracy. All methods represent the region be- 
hind the contact discontinuity and the expansion wave 
equally well and no major differences exist between the 
reversed or natural ordering of the eigenvalues in the 
Osher scheme. It should be borne in mind that due 
to the fact that the solution of the Riemann problem 
is self similar only first order accuracy is possible and 
the shock tube tests merely show that solutions with 
higher order methods can be obtained without numer- 
ical oscillations. Test of all methods on free convection 
problems, however, showed that all methods reached 
the proper order of accuracy. Higher order methods, 
especially ENO-CK, however, still give a better solu- 
tion than first order methods for flows with disconti- 
nuities. Especially the smearing in the region between 
the shock and contact discontinuity reduces. 

The second test, case B, is different from case A in that 
it has a sonic point in the expansion wave. The first or- 
der solution therefore has an O(Az) expansion wave at 
the sonic point, but does not need an entropy fix as does 
the Roe method. This discontinuity is stronger for the 
natural order than for the reversed order of eigenvalues 
in the Osher scheme. The disturbance is about twice 
as large for the natural order and although greatly re- 
duced for higher order methods it consistently slightly 
reduces the accuracy in this area, when using the nat- 
ural order. This is a general conclusion for all cases 
with a sonic point. Both ENO-CK and ENO-FL work 
reasonably well, with ENO-CK, with reversed order of 
eigenvalues, superior, as can be seen in Fig. 4-5. The 
flux-ENO method experiences a small jump around the 
sonic point due to the limited differentiability of the Os- 
her flux at this point. This problem would require ad- 
ditional attention when using the flux-ENO method to 
obtain uniform high order solutions, but for the com- 
parison of different ENO methods it is not relevant. 
The ENO-RK method experienced serious oscillations 
at the beginning and end of the expansion wave, Fig. 6, 
and was unstable for fifth order accuracy, when using 
the natural ordering. Both natural and reversed order 
became oscillatory, but the natural order turned out to 
be more sensitive in general. In cases of instability the 
differences, however, were not large enough to prefer 
one method for the other. The other two methods also 
experience small problems at the bottom of the expan- 
sion wave, but much less severe. A clear explanation for 
this phenomenon still is lacking and will require further 
research. 

The third problem, case C, has as main feature a strong 
contact discontinuity. The different 


Order 

ENO-FL 

ENO-CK 

ENO-RK 

2 

1 . 

1.13 

2.15 

3 

1.70 

1.73 

3.68 

4 

1.91 

2.07 

4.15 

5 

2.13 

3.12 

4.62 


Table 2. Comparison CPU time of different 
ENO schemes 

schemes behave similar in the region with the expan- 
sion wave as in case B, with the reversed order slightly 
better around the sonic point. All methods experience 
a small undershoot in the region between the contact 
discontinuity and the expansion wave, with ENO-CK, 
Fig. 7, the most accurate and ENO-FL has some high 
frequency oscillations, see Fig. 8. The smearing of 
the contact discontinuity is approximately equal for all 
methods, but slightly reduces with higher order accu- 
racy. Contrary to case B, ENO-RK now does not ex- 
perience instability problems. 

The final test, case E, which consists of two strong 
shocks moving in opposite directions is the most dif- 
ficult case. The first order Osher scheme with natural 
ordering of eigenvalues is not monotone and has a sig- 
nificant overshoot in density, pressure and velocity at 
the first point behind the left moving shock, Fig. 9. 
The Osher scheme with reversed ordering is monotone. 
All ENO schemes are already oscillatory for second or- 
der accuracy. The second order ENO-CK scheme, with 
reversed ordering of eigenvalues experiences some small 
oscillations in the pressure behind the shock, Fig. 11, 
but the density still is nearly monotone, Fig. 10. The 
reason for the oscillations are the same as those dis- 
cussed for case A. In the initial stages there are not 
enough points between the two shocks to build a non- 
oscillatory reconstruction. This phenomenon is, how- 
ever, stronger for case E than for case A, because this 
problem only occurred in case A for the density which 
had a region separated by a shock and contact discon- 
tinuity, but not for the other variables. When a region 
is separated by two shocks also the pressure and veloc- 
ity are discontinuous. Contrary to the expectations the 
second order flux-ENO method, which was oscillatory 
for case A now is nearly monotone. All ENO schemes 
are oscillatory for third and higher order, as can be seen 
in Fig. 12-13. 

The choice between the different ENO methods is not 
only determined by its accuracy but also by its numer- 
ical efficiency. Table 2. shows the relative CPU time 
for the different methods. The flux-ENO method is 
the least expensive, but the ENO-CK scheme is nearly 
equal up to 4-th order. The difference at 5-th order 


is caused by the fact that the other two schemes only 
use the third order accurate TVD Runge-Kutta scheme, 
while ENO-CK also is fifth order accurate in time. Ta- 
ble 2 also shows that the ENO-RK scheme is signifi- 
cantly more expensive and requires approximately twice 
as much CPU time. 

5. Conclusions 

In all cases the ENO-CK method, with reversed order- 
ing of eigenvalues in the Osher scheme, was superior 
or performed equally well as the other ENO schemes 
and is the most robust. The ordering of eigenvalues in 
the Osher scheme is not extremely critical. The natu- 
ral ordering, however, has a larger jump around a sonic 
point, which slightly reduces the accuracy of higher or- 
der approximations. The natural ordering also expe- 
riences a slightly larger level of numerical oscillations 
than the reversed ordering. In cases were the first order 
scheme with natural ordering is not monotone it is not 
possible to build higher order non-oscillatory schemes. 
This stresses the importance of a good approximate 
Riemann solver. 

The ENO method, with TVD Runge-Kutta time inte- 
gration, is more dissipative and less robust than the 
ENO method with the Cauchy-Kowalewski procedure 
for time integration. Although it requires a signifi- 
cant effort to derive all the higher order derivatives 
for a multi-dimensional problem, it certainly pays off 
when one considers the fact that the TVD Runge-Kutta 
method requires multiple reconstructions, which are the 
most expensive part of the algorithm, and significant 
larger storage. The flux-ENO method, although hav- 
ing the benefit of being the most easy to program, and 
fairly straightforward to extend to multiple dimensions 
does not have the robustness of the primitive variable 
reconstruction. 
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Figure 7. Case C, density at t — .5, 2nd (dots) and Figure 10. Case E, density at i — .6, 2nd order ENO- 

5th order (triangles) ENO-CK. CK (dots) and ENO-FL (triangles). 




Figure 8. Case C, density at t = .5, 2nd (dots) and Figure 11. Case E, pressure at i =. .6, 2nd order 

5th order (triangles) ENO-FL. ENO-CK (dots) and ENO-FL (triangles). 



Figure 0. Case E, density at t — .6, 1st order Osher 
scheme, natural (triangles) and reversed order (dots). 
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Figure 12. Case E, density at i — .6, 3nd order ENO- 
CK (dots) and ENO-FL (triangles). 


0.1 



Figure 13. Case E, pressure at t — . 
ENO-CK (dots) and ENO-FL (triangles). 
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